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Abstract 


The Crau Plain grasshopper, Prionotropis rhodanica Uvarov, 1923 (Or- 
thoptera: Pamphagidae: Thrinchinae), is a rare grasshopper species en- 
demic to the Crau Plain, a steppic habitat in France with unique floristic 
and faunistic communities. During recent decades, the area covered by 
these steppic grasslands has been highly reduced and fragmented due to 
the development of irrigation-based agriculture, roads, as well as industrial 
and military complexes. The restricted distribution, low population den- 
sity and poor dispersal ability of P. rhodanica, combined with the destruc- 
tion of its habitat, has led to the classification of this species as critically 
endangered in the IUCN Red List of Threatened Species. Decreases in habi- 
tat quality due to intensive grazing in the remnant grassland patches con- 
stitute an additional threat for P. rhodanica that can impact population dy- 
namics at a relatively small-scale. In this work, we focused on a small area 
of about 3 km? occupied by one of the largest subpopulations observed 
in 2000-2001. We conducted a single-time snapshot intensive survey of 
grasshopper density and genetic variation at 11 microsatellite markers. We 
used a recent method, MAPI, to visualize the spatial genetic structure as a 
continuous surface and to determine, with the simultaneous use of spatial 
cross-correlograms, whether the normalized difference vegetation index, 
which informs on the balance between vegetation productivity and graz- 
ing intensity, can explain grasshopper population structure at such a fine 
scale. We found that both population density and gene flow were strongly 
and positively correlated to habitat quality (higher productivity of grass- 
lands and/or lower sheep grazing). The spatial scales of interaction be- 
tween these variables were estimated to be highly similar, in the range of 
812-880 meters. This result suggests that P. rhodanica is very sensitive to 
the quality of the grasslands it inhabits. 
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Introduction 


The Crau Plain (Bouches du Rhéne, France; Fig. 1A-B), the 
ancient delta of the Durance River, is the last arid steppe in France 
with a unique fauna and flora (Cheylan 1975, Devaux et al. 1983, 
Wolff et al. 2001, 2002, Romermann et al. 2005, Tatin et al. 2013). 
The characteristic vegetation of this steppe is called Coussoul and 
is mainly composed of Brachypodium retusum and Thymus vulgaris 
in association with Asphodelus fistulosus and Stipa capillata. This 
unique ecosystem has been maintained over 3,000-4,000 years 
by extensive sheep grazing (Buisson and Dutoit 2006, Tatin et 
al. 2013). Since the sixteenth century, and more intensely within 
recent decades, the plain has been highly modified by intensive 
irrigation-based agriculture (Devaux et al. 1983). The contraction 
of the Coussoul from ~600 km? to less than 100 km? (Wolff et al. 
2002) and its fragmentation by agricultural fields (orchards and 
hay meadows), irrigation channels and other artificial structures 
(roads, industrial and military complexes) threaten the natural 
habitat and its fauna. Since 2001, 7,400 ha of the Crau Plain have 
been protected as a national nature reserve and since 2010 a man- 
agement plan has been implemented with conservation actions 
for its remarkable species. Among them, P. rhodanica (Fig. 2) is a 
large (between 30 and 50 mm) grasshopper species endemic to 
the Crau Plain and protected in France (Anonyme 1979). This spe- 
cies has been classified as Critically Endangered in the IUCN Red 
List of Threatened Species (Hochkirch and Tatin 2016) as well as 
on the European Red List (Hochkirch et al. 2016) and the national 
Red List of France (Sardet and Defaut 2004). P. rhodanica is only 
known from the Crau Plain, though surrounding areas have been 
intensely surveyed. The species is generally considered as rare, ex- 
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Fig. 1. Map of France and location of A. the plain of Crau (Bouches-du-Rh6ne, France), B. the sampling site located between the sheep- 
folds ‘La Grosse du Levant’ and ‘La Grosse du Centre’, and C. the circles surveyed to detect and sample P. rhodanica. AO0-A65: names 
of circles of the grid A; B04-B51: names of circles of the grid B; C1: additional circle of 100 m diameter (see Material and methods). 


JOURNAL OF ORTHOPTERA RESEARCH 2018, 27(1) 


S. PIRY, K. BERTHIER, R. STREIFE, $. CROS-ARTEIL, A. FOUCART, L. TATIN, L. BRODER, A. HOCHKIRCH AND M.-P. CHAPUIS 63 


Fig. 2. Adult of Prionotropis rhodanica. 


cept for some atypical years, and observations over the last fifteen 
years revealed an extreme decline in population densities along 
with local extinctions (Foucart and Lecoq 1996, Foucart et al. 
1999, Hochkirch et al. 2014). In 2014, a conservation strategy for 
the species was developed under the guidance of the Species Con- 
servation Planning Sub-Committee (SCPSC) and the Invertebrate 
Conservation Sub-Committee (ICSC) of the International Union 
for Conservation of Nature (IUCN) (Hochkirch et al. 2014). One 
main aim of the strategic conservation plan was to identify the 
reasons for the population decline. 

Orthoptera are known to be highly sensitive to landscape al- 
terations, such as those associated with farming, in terms of genet- 
ic diversity and structure (Ortego et al. 2012, Gauffre et al. 2015, 
Ortego et al. 2015) as well as population decline or extinction risk 
(Barker 2004, Reinhardt et al. 2005). However, demographic re- 
sponses and evolutionary trajectories differ among species, and 
some species are more susceptible than others regarding the nega- 
tive effects of habitat perturbations (Ortego et al. 2015). Thus, 
conservation practices in protected areas require detailed ecologi- 
cal and evolutionary information to target species that require par- 
ticular attention and guide their management (Ortego et al. 2015). 
A capture-mark-recapture study conducted on P. rhodanica at a fine 
scale (8,000 m?) over a two-month period (May-June) in 2000 
showed that individual movements occur over very short distanc- 
es (mean distance between two captures < 20 m) and suggested 
that adults are more mobile than nymphs (Berthier 2000, see also 
Besnard et al. 2007). These very limited dispersal abilities of P. 
rhodanica were expected as adults of both sexes are brachypterous 


and thus incapable of flight (Fig. 2; Foucart 1995, Foucart et al. 
1999). Furthermore, studies of the genetic structure of P. rhodanica 
conducted on a small number of sample sites (2-5) and micros- 
atellite loci (5) revealed that subpopulations are strongly differ- 
entiated even when located a few hundred meters apart (Berthier 
2000, Streiff et al. 2005). These results suggest that P. rhodanica 
subpopulations are isolated with dispersal events among them be- 
ing rare, confirming that biological and ecological characteristics 
of the species result in low demographic and genetic connectiv- 
ity between isolated subpopulations. Isolation effects depend on 
the sizes of habitat elements, the distances between them and the 
quality of the surrounding matrix (Prevedello and Vieira 2010). 
Their quantification would thus be important to assess the species 
vulnerability to the reduction and fragmentation of its habitat and 
would require an intensive sampling of remnant subpopulations 
at a large scale (e.g. the complete geographic range of the species). 

An additional threat for P. rhodanica might result from live- 
stock grazing pressure (Hochkirch et al. 2014). While the available 
surface of Coussoul has continuously been decreasing for centu- 
ties, flock sizes fluctuate from one year to the next and can reach 
1,600 sheep on average (Wolff et al. 2013). Moreover, the size of 
flocks can strongly increase when sheep are gathered to prepare for 
the transhumance in June-July, coinciding with the breeding peri- 
od of P. rhodanica (Foucart and Lecog 1996). For example, during 
the capture-mark-recapture study conducted in 2000, the size of 
the flock grazing the 3 km? land plot within which our study site 
of 8,000 m? was localized increased from 250 to 1,200 individuals 
in two months (Berthier 2000). Generally, such successive peri- 
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Fig. 3. Illustrations of the Coussoul habitat of Prionotropis rhodanica with A. sheep accompanied by cattle egrets and B. researchers sur- 


veying the species using a circle-based method. 


ods of intensive grazing can be detrimental because they reduce 
primary productivity, impede plant growth and alter vegetation 
cover that provides the grasshopper with food and shelter. In P. 
rhodanica, the highest observed abundances were associated with 
a vegetation cover above 70% (Tatin et al. 2013). In addition to 
competition for resources, negative impacts of intensive grazing 
include trampling and predation by birds associated with sheep 
flocks (Hochkirch et al. 2014). In particular, the cattle egret, Bubul- 
cus ibis L., often accompanies sheep flocks (Fig. 3A) and catches 
insect and small vertebrate prey disturbed by sheep (personal 
observation). Since pasture boundaries delineate plots less than 
550 ha (min of 30 ha and mean of 168 ha), grazing pressure and 
consequently grasshopper habitat quality may vary at a relatively 
small scale in the Crau Plain, even within apparently continuous 
and favorable areas of Coussoul. Information on fine-scale spatial 
genetic structure would allow an assessment of the role of habitat 
quality, rather than that of unsuitable elements that isolate habitat 
patches, on contemporary dispersal patterns. However, a sampling 
scheme at a fine scale has never been attempted in P. rhodanica. 


MAPI (Mapping Averaged Pairwise Information; Piry et al. 
2016) is a new method for spatial visualization of population 
structure and investigating landscape effects. MAPI is essentially a 
smoothing procedure for spatial genetic networks, useful for large 
sample sets for which usual representations (i.e. nodes and edges) 
are often unreadable. MAPI can be applied to genetic differentia- 
tion measures computed from neutral markers (e.g. microsatel- 
lites) to detect areas of high genetic continuity and/or discontinuity 
that reflect areas where gene flow is the highest and/or the lowest, 
respectively. MAPI allows visualizing the spatial genetic structure as 
a geographic map that provides information on the average intensi- 
ty of the genetic relationships, and can be exported simultaneously 
with landscape layers and density data for further statistical analy- 
ses. This exploratory approach may thus provide information on 
which environmental variables are potential candidates to explain 
observed genetic patterns. A central feature of MAPI is its low sen- 
sitivity to confounding effects resulting from isolation-by-distance 
(IBD). This is especially true in an ideal situation, with perfect regu- 
lar sampling and no edge effects (Piry et al. 2016). However, when 
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sampling is highly irregular, the insensitivity to IBD still holds rela- 
tively to clustering methods (Guillot and Santos 2009, Bradburd et 
al. 2013, Piry et al. 2016). This is a critical issue when assessing the 
impact of landscape heterogeneity on population genetic structure, 
which also depends on geographic distance (Bradburd et al. 2013). 
This may be of high relevance to P. rhodanica since significant IBD 
patterns have been reported at fine spatial scale (2-10 km) in other 
orthopteran species that are flightless (Lange et al. 2010, Keller et al. 
2013, Gauffre et al. 2015), or considered to be sedentary (Ortego et 
al. 2011, Blanchet et al. 2012). 

In this study, we focused on a small area of 2.87 km? within 
which direct gene flow was likely to occur while potential effects of 
the physical landscape (i.e. barriers to dispersal, presence of inhos- 
pitable patches) were expected to be minimal. Our target study site 
was a Suitable habitat for P. rhodanica but with potential variation 
in quality in particular due to sheep grazing regime. It was located 
in a large patch of Coussoul where the largest subpopulation of 
P. rhodanica was observed at the end of the nineties (Foucart et al. 
1999). In spring 2001, we carried out a single-time snapshot inten- 
sive survey of grasshopper density and genetic variation at eleven 
microsatellite loci. The sampled subpopulation was composed of 
266 individuals of a single cohort since P. rhodanica is univoltine 
with discrete generations (eggs hatch in early April and all adults 
die by mid-July; Foucart and Lecoqg 1996). We investigated con- 
temporary patterns of dispersal by first testing whether there is 
non-random spatial structure due to relatedness among individu- 
als at the scale of the sampling unit (i.e. < 50 m) and/or isolation- 
by-distance at the scale of the study area. We then used MAPI to 
produce a continuous variation surface of the genetic relationships 
between individuals. Finally, we used spatial cross-correlograms to 
quantify the spatial scales of interactions between this genetic pat- 
tern and variations in grasshopper density and habitat quality. To 
quantify small-scale variation in habitat quality in our study site, 
we used the normalized difference vegetation index (NDVI) as a 
proxy of the balance between vegetation productivity and vegeta- 
tion degradation partly due to grazing. Remotely sensed vegeta- 
tion indices are commonly used to monitor changes in vegetation 
and land cover, and their correlation with grazing intensity has 
been demonstrated in other parts of the world (e.g. Raynolds et al. 
2015), including in semi-arid environments (Blanco et al. 2008, 
Karnieli et al. 2013). Here, the validity of the relationship between 
vegetation productivity, NDVI and grazing effect was evaluated in 
the specific area of the Crau Plain using data from a recent experi- 
ment (i.e. conducted in 2015). 


Material and methods 


Sampling and grasshopper densities. —P. rhodanica is particularly dif- 
ficult to detect in the field (less than 10% detection probability) 
due to its low mobility, cryptic coloration and crypsis behavior 
(Streiff et al. 2005, Besnard et al. 2007, Tatin et al. 2013). There- 
fore, we determined a sampling area of 2.87 km? to focus on the 
zone where densities were relatively large in 2001 (Fig. 1B-C). 
Within the study area, circles with 50 m diameter were regularly 
distributed to create two completely overlapping grids (A and B), 
so that circles from grid A alternate with those from grid B eve- 
ry 160 meters (Fig. 1C). All circle centers were positioned using 
Global Positioning System (GPS). Grid A was entirely surveyed be- 
tween the 21% and 29" of June 2001 to ensure a complete coverage 
of the study area. In grid B, only the circles located in zones where 
the species was detected in grid A were surveyed between the 27th 
to the 31st of June 2001. In addition to the 50 m circles of grid A 


and B, a circle with 100 m diameter located within the zone with 
the highest P. rhodanica density was also surveyed between the 21* 
and 29" of May 2001. 

All circles were divided in 8 adjacent sectors (16 for the 100 m 
diameter circle), physically delimited in the field, and each sector 
was explored by three persons walking abreast and passing at least 
twice on a same path (Fig. 3B). For each insect observed, polar coor- 
dinates to the center were recorded, and then converted in latitude 
and longitude. Sex and developmental stage (nymph/adult) were 
recorded. A non-destructive sampling method was used for further 
genetic analyses (a small piece of elytra or tarsus was collected and 
the insect released at its capture position). From the 82 surveyed 
circles, a total of 266 individuals (213 nymphs, 53 adults, 132 fe- 
males and 134 males) were sampled. More than 50% of the circles 
(i.e. 44) did not contain grasshoppers and densities varied from 5.1 
to 133.3 individuals per hectare within the 38 remaining circles 
(Suppl. material 1: Table $1). Accordingly, an analysis of the spatial 
distribution using Ripley’s K function revealed a patchy distribution 
of the individuals within the study area (Suppl. material 1: Fig. §1). 


Microsatellite genotyping and basic genetic indices. -Genomic DNA was 
extracted following the CTAB protocol (Doyle and Doyle 1987). 
DNA concentration was determined using a Nano-drop spectropho- 
tometer (NanoDrop Technologies, Wilmington, USA) and extracts 
were normalized to a concentration of ca. 5 ng/L. Eleven microsat- 
ellite loci were used, ten of which were previously described in Streiff 
et al. (2002) and Streiff et al. (2005) and one additional locus was 
developed for this study (see primer sequences in Suppl. material 1: 
Table S2). We genotyped the 266 individuals at the 11 microsatellites 
using an ABI PRISM 310 DNA sequencer (Applied Biosystems) us- 
ing conditions described in Suppl. material 1: Table S2. Alleles were 
scored using GeneMapper 4.0 (Applied Biosystems). 

The level of polymorphism and allelic distribution were es- 
timated with GENEPOP v.4 (Rousset 2008). We first tested for 
linkage disequilibrium between each pair of loci and within each 
circle by using G-exact tests (keeping the 22 circles out of 38 
non-empty circles with sampling size > 3 individuals). Over the 
sampling site, and for each locus, we tested for Hardy-Weinberg 
equilibrium using G-exact tests and estimated the number of al- 
leles, observed and expected heterozygosities of Nei (1987) and 
Wright's F-statistic F., according to Weir and Cockerham (1984). 
We also tested for genotypic differentiation among circles by using 
the exact G-test of heterogeneity of genotypic frequencies (keeping 
the 22 circles with sampling size > 3 individuals). 


Kinship and isolation-by-distance patterns.—In order to analyze the 
relatedness structure at the circle scale, we calculated the kinship 
coefficient of Loiselle et al. (1995) between pairs of individu- 
als belonging to a same circle, using SPAGeDI v.1.4 (Hardy and 
Vekemans 2002). We then used Spearman rank correlation to 
test whether there is an association between grasshopper density 
measures and mean values of kinship coefficient within circles. 
We assessed whether dispersal was restricted with distance, 
using GENEPOP v.4 (Rousset 2008). We computed for each pair 
of sampled individuals, and over all loci, the statistics 4 (Rousset 
2000), which is considered as an unbiased estimator of genetic 
distance (Watts et al. 2007). We used a Mantel test between the 
matrices of genetic distances (4) and the logarithm of geographi- 
cal Euclidean distances to test for a positive linear relationship as 
expected under an IBD model. We excluded pairwise comparisons 
between individuals separated by less than 75 m to exclude intra- 
circle measures. This would limit the impact of siblings on the es- 
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timation of the relationship between genetic and geographic dis- 
tances. Confidence interval and significance of regression slope was 
calculated by bootstrapping over loci using 10,000 permutations. 
We also performed the same IBD analysis for each sex separately. 


NDVI computation.—As suggested by Gan et al. (2014), we com- 
puted rescaled NDVI values by building fine-scale NDVI raster 
from Landsat images and calibrating the results using a coarse- 
resolution MODIS raster. The LANDSAT 7/8 ETM+ raster (30 m 
resolution) and NDVI MODIS raster, all courtesy of the U.S. Geo- 
logical Survey, were downloaded from https://earthexplorer.usgs. 
gov. The MODIS NDVI is a 16-day composite of MODIS data at 
a spatial resolution of 250 m. First, NDVI values were computed 
at the spatial resolution of the Landsat raster (30 m) using an up- 
dated version of the python script provided by Dr. J. Degener and 
available from https://www.uni-goettingen.de/en/524379.html. 
Second, the Landsat NDVI raster and the MODIS NDVI raster were 
clipped over the area of interest (projection UTM 31 / WGS84). 
Third, the MODIS NDVI raster was aligned on the Landsat NDVI 
raster using the nearest value. Both rasters were then loaded in 
the statistical software R (R Core Team 2015) using the packages 
“tgdal” and “raster” and, finally, the pixel values from both rasters 
were bound into a single dataframe. Pixels with missing MODIS 
values (-3,000) were discarded from the dataset. MODIS NDVI 
values were divided by 10,000 to be within the range of the Land- 
sat NDVI values. Landsat NDVI values were then normalized by 
performing a pixel-to-pixel regression of the MODIS NDVI values 
against the Landsat NDVI values (see Suppl. material 1: Figs $2, $3 
for an illustration). 


Relationship between NDVI, habitat quality and grazing.—We ana- 
lyzed the relationship between rescaled NDVI, vegetation produc- 
tivity and grazing using data from a recent experiment. In 2015- 
2016, an enclosure of 8.56 ha, located in the same area as our 
study site (Suppl. material 1: Fig. $4), was fenced between April 
and August to exclude sheep. One hundred plots (of 30 cm diam- 
eter each) were randomly selected to measure vegetation produc- 
tivity indices within (50 plots) and outside (50 plots) this enclo- 
sure. Vegetation productivity indices included: vegetation height, 
coverage in forbs, coverage in dry vegetation and coverage in bare 
ground (i.e. no vegetation). We tested whether our four produc- 
tivity indices reflected the impact of sheep grazing on vegetation 
using Wilcoxon tests to compare the measures between plots lo- 
cated inside and outside the enclosure. Rescaled NDVI values were 
computed from Landsat and MODIS images captured in May 2014 
(before fencing) and May 2016 (after fencing). Each vegetation 
plot received the rescaled NDVI value of the 30 m pixel to which 
it belonged. We analyzed the relationship between 2016-rescaled 
NDVI values (after fencing) and the four vegetation productivity 
indices using Spearman rank correlations. Finally, we used Wilcox- 
on tests to verify that after fencing (2016) the rescaled NDVI values 
were, on average, significantly lower for the plots located inside 
the enclosure than outside, while no differences were expected for 
2014-rescaled NDVI values (before fencing). In other words we 
tested that our rescaled NDVI was an appropriate proxy to assess 
grazing impact on grasshopper habitat quality. Statistical analyses 
were done using R (R Core Team 2015). 


Spatial association between grasshopper density, habitat quality and 
genetic variation.—We used the MAPI program (Piry et al. 2016), 
freely available at https://www1.montpellier.inra.fr/CBGP/soft- 
ware/MAPI/. We produced a geographical map of the spatial vari- 


ation of the mean genetic differentiation between individuals and 
then interpolated, within the cells of the map grid, grasshopper 
density and habitat quality (rescaled NDVI) to test for correla- 
tions between variables. The successive three steps of this analysis 
are detailed below. 


Map of genetic differentiation.—We attributed the values of the sta- 
tistics @ of Rousset (2000) to ellipses that spatially materialize 
the connections between the pairs of individuals. A grid of hex- 
agonal cells with a half-width of 20 m was superimposed on the 
study area (cell’s area = 1,040 m/’; see ‘Grid of cells’ section in Piry 
et al. (2016) for details). Each cell received the weighted arith- 
metic mean of the ellipses intercepting its geographical extent. 
This mean was weighted using the inverse of the ellipse areas 
to limit long distance effects. MAPI values for the eccentricity of 
ellipses that controls the smoothing intensity was set by default 
(i.e. eccentricity = 0.975). However, additional analyses were run 
with eccentricity values from 0.800 to 0.999 to assess robustness 
of MAPI results to the shape of ellipses (i.e. from inflated to nar- 
row). The radius of the error circle that controls for uncertainty 
on sample polar coordinates (error_circle_radius) was set to 2 m. 
An interesting feature of MAPI is the optional parameters that 
limit the analysis to a given range of between-sample distances. 
We used a minimum distance between samples (min_distance) 
of 75 m in order to exclude pairwise comparisons between indi- 
viduals from a same circle. Details on parameters can be found 
in Table 1 and Suppl. material 1 from Piry et al. (2016) and in 
pp. 29-31 from the MAPI manual. Finally, to ensure that MAPI 
results were not biased by the patchy distribution of individuals 
(higher number of individuals sampled in high density areas), 
we ran the analysis on 100 independent resampled datasets with 
a maximum of three individuals from each circle and checked for 
consistency of results. 


Density and NDVI spatial interpolation.—Grasshopper densities 
measured within circles and rescaled NDVI values from Landsat 
(captured the 18 of May 2001: LEO7_L1TP_196030_20010518_2 
0170205_01_T1.tar.gz) and MODIS (captured the first 16 days of 
May 2001: MOD13Q1.A2001129.h18v04.006.2015142065654. 
hdf) images were first interpolated into a raster (5 m resolution) 
using the function “using v.surf.rst” of the GRASS software (GRASS 
Development Team 2012) with parameter values set as follows: 
maximum number of points in a segment (segmax) = 120, mini- 
mum number of points for approximation in a segment (npmin) 
= 60 and maximum distance between points on isoline (dmax) = 
25. Values interpolated below 0 were reset to 0 using the GRASS 
function “rmapcalc”. Second, raster pixels were clipped to MAPI 
grid cells and interpolated density values from pixels belonging 
to the same MAPI cell were averaged using the ST_SummaryStats 
function of the PostgreSQL 9.3 with the PostGIS 2.1.2 extension 
(1996-2016, The PostgreSQL Global Development Group: http:// 
www.postgresql.org/ — http://postgis.net/). 


Spatial scales of association between variables.—Spatial cross-correlo- 
grams allow investigation of how two variables co-vary with geo- 
graphic distance. We used the non-parametric spline-correlogram 
approach implemented in the R package “ncf” (Bjgrnstad and Fal- 
ck 2001) to analyze the spatial scale of association (Sp) between 
interpolated values of: 1) grasshopper density and rescaled NDVI, 
2) grasshopper density and mean genetic differentiation between 
individuals (i.e. MAPI cell values) and, 3) rescaled NDVI and mean 
genetic differentiation between individuals (i.e. MAPI cell values). 
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Confidence envelopes at 95% (CE,,,,) for the estimated correlo- 
grams were calculated using bootstrapping (500 replicates). We 
also computed Spearman rank correlation coefficients (Rho) be- 
tween all pairs of variables. 


Results and discussion 


Thirty-four of the 1,210 tests for linkage disequilibrium be- 
tween the 11 loci were significant after false discovery rate correc- 
tion (Storey and Tibshirani 2003). Because these pairs of loci were 
never significant in more than two circles, all microsatellite loci 
were considered unlinked. We found significant deviations from 
Hardy-Weinberg equilibrium at the global scale for most loci and 
the F,, value averaged to 0.072 across loci (Table 1). Heterozygote 
deficits can be partly related to presence of null alleles for the few 
loci that show highest values (i.e. Phr7178 and Phr756 with a F., 
> 0.3). Prevalence of null alleles has commonly been reported 
in studies documenting microsatellite variation in orthopteran 
populations (Zhang et al. 2003, Chapuis et al. 2005, Hamill et 
al. 2006, Berthier et al. 2008, Chapuis et al. 2008). However, the 
heterozygote deficit concerns most loci, which suggests that it pri- 
marily results from spatial structure (i.e. Wahlund effect) rather 
than presence of null alleles. As a matter of fact, global genetic 
differentiation between circles was highly significant (P < 0.0001) 
and all pairs of circles were significantly differentiated (P < 0.001). 
Moreover, most of the Loiselle kinship coefficient values averaged 
within circles were positive. Thus, pairs of individuals within cir- 
cles were more related than expected under random distribution 
of genotypes. This pattern at a scale of 50 m confirmed the limited 
dispersal capacities of this flightless grasshopper species, at least at 
the nymphal stage. 

The levels of genetic diversity were high, with a mean expected 
heterozygosity of 0.813 and an average of 16.18 alleles for our 
sample size of 266 individuals (Table 1). They are similar to the 
levels of genetic diversity observed by Berthier (2000) and Streiff 
et al. (2005) within subpopulations of the same species. The 
(apparent) lack of impact of habitat modifications on the level 
of genetic diversity may be explained by their recentness in the 
Crau Plain (Streiff et al. 2005). Changes in genetic diversity after 
a perturbation are related to the total effective population size, 
and show slow dynamics with non-equilibrium states over large 
temporal scales. This result may also account for the limited dis- 
persal in the fragmented population of P. rhodanica, which allows 
for strong differentiation between subpopulations and thereby a 
high level of total genetic diversity (Nei 1973). Beside population 
demography and history, this outcome may also result from the 
high microsatellite loci mutation rate, which has been confirmed 
for the Orthoptera (Chapuis et al. 2012). Overall, P. rhodanica did 
not seem to be affected by a low level of genetic diversity, which 
is widely recognized as a major impediment for the adaptation of 
a population to environmental changes. However, this does not 
preclude that the species may have a low total population size, 
which would make it vulnerable to demographic stochasticity that 
can lead to local extinctions (Frankham 2005). Indeed, during a 
standardized survey conducted in 2012-2013, P. rhodanica was not 
detected in our study site anymore and this subpopulation is now 
assumed to be extinct (Hochkirch et al. 2014). 

We detected a significant negative relationship between the 
Loiselle kinship coefficient and density within circles (Rho = 
-0.59; p-value = 0.0012; Fig. 4A), ie. the lower the grasshopper 
density, the higher the genetic relatedness. Since we did not find 
any evidence for lower levels of genetic diversity within lower 


Table 1. Genetic diversity indices for each of the 11 loci from P. 
rhodanica. Number of alleles (Na), Wright's F-statistic F,, and ob- 
served (Ho) and expected (He) heterozygosities were averaged 
over the 266 individual samples. 


Locus Na Be Ho He es 
P value 
Phric7 23 0.034 0.865 0.895 < 0.0001 
Phr228 9 0.028 0.545 0.561 OA'770 
Phr2C3 10 0.015 0.816 0.828 < 0.0001 
Phr2T 12 -0.021 0.774 0.758 0.0279 
Phr3B3 23 0.031 0.857 0.884 0.0238 
Phr4A10 34 0.005 0.936 0.941 < 0.0001 
Phr4G1 keg 0.142 0.774 0.901 < 0.0001 
Phr7178 5 0.316 0.432 0.632 < 0.0001 
Phr756 0.306 0.510 0.734 < 0.0001 
Phr880 19 -0.037 0.940 0.906 OCSEs i531) 
Phr4H3b i 0.059 0.853 0.907 < 0.0001 
Over all loci 16.18 0.072 0.755 0.813 < 0.0001 


density circles, we can exclude the effects of genetic drift under 
random mating as cause (Rho = +0.14; p-value = 0.506; Fig. 4B). 
Moreover, negative F,, values were associated with the lowest grass- 
hopper density conditions and thereby with the highest genetic 
relatedness values (Rho = +0.62; p-value = 0.0019; Fig. 4C). Theo- 
retical studies demonstrate that excess heterozygosity is expected 
when allelic frequencies differ in fathers and mothers (e.g. Rousset 
and Raymond 1995). For example, in social species negative val- 
ues of F,, may be indicative of complicated breeding tactics rather 
than classical outbreeding (Van Staaden 1995) as exemplified in 
some populations of elephants, bats and finches (Tarr et al. 2000, 
Nyakaana et al. 2001, Storz et al. 2001). Since grasshoppers are 
often polyandrous, our result can hence be explained by the fact 
that most of the nymphs sampled from the lowest density circles 
may have emerged from a single egg-pod, where half-siblings are 
from a same mother but several unrelated fathers. 

We detected a significant positive linear relationship between 
the differentiation coefficient (4) and the logarithm of geographi- 
cal distance at a scale < 2,500 m (P < 0.0001; Fig. 5). The value 
of the regression slope suggests that dispersal decreased strongly 
with geographical distance in this species (estimate [95% confi- 
dence interval] = 0.017 [0.008-0.033]). Significant IBD patterns 
were still found when analyzing males and females separately and 
overlapping confidence intervals for slope estimates did not sup- 
port sex-biased dispersal, at least at the nymphal stage (Fig. S5 in 
the Suppl. material 1). Overall, this result indicates that dispersal 
is seriously restricted in space at the sampled scale. This dispersal 
pattern is likely to limit the ability of the species to colonize new 
areas and can ultimately reduce the long-term persistence of the 
isolated populations. This is in agreement with local extinctions 
observed since 1996 in areas that have never been recolonized 
since (Foucart and Lecoq 1996, Foucart et al. 1999). 

The map of the interpolated densities visually confirmed the 
result of the Ripley’s K statistics with a clear occurrence of two main 
high density nuclei in the northern half of the study site, while the 
density was very low in the southern half (Fig. 6A). An equivalent 
pattern was observed from the rescaled NDVI map, with highest 
values occurring in the northern half of the study site (Fig. 6B). 
The MAPI analysis also revealed a strong spatial structure, with the 
lowest levels of genetic differentiation observed between individu- 
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als located within the high density nuclei (Fig. 6C). Individuals lo- 
cated in lower density areas (extreme southern and northern parts 
of the study area) were not only differentiated from the density 
nuclei but also exhibited a relatively high level of differentiation 
across short distances. It is worth noting here that these results are 
expected to be robust to the observed pattern of strong IBD, even 


Fig. 4. Relationships between intra-circle grasshopper density and A. Loiselle kinship coefficient, B. 


B43" 


expected heterozygosity and C. F... 


under our irregular sampling data (Piry et al. 2016). In addition, 
the optional parameter of minimum distance between samples in 
MAPI allowed us to eliminate the potential impact of the inferred 
presence of siblings within circles of low grasshopper density on 
the assessment of the surface of genetic variation. MAPI results 
from resampled datasets (of a maximum of 3 individuals within 
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Fig. 5. Linear regression between genetic distance 4 and geographical distances computed between pairs of individuals. Variation in 
point density is represented by colors, from blue (low density) to red (high density). 


circles) were highly consistent with those obtained from the com- 
plete dataset (Suppl. material 1: Fig. S8). 

Finally, Spearman coefficients showed that grasshopper den- 
sity and rescaled NDVI values were positively correlated (Rho = 
0.51; p-value < 2.2e'°) but both negatively correlated to MAPI 
cell values (Rho = -0.42 and -0.34, respectively; p-value < 2.2e 
‘6 for both). The three cross-correlograms showed that the spa- 
tial scales of association between variables were highly similar 
and quite small (Fig. 7): rescaled NDVI and density values, Sp = 
880m (CE,.,, = [855-902]); density and MAPI cell values, Sp = 859 
(CE,,,, = [680-927]); rescaled NDVI and MAPI cell values, Sp = 
812m (CE,,,, = [781-850]). The strong association between values 
of grasshopper density, rescaled NDVI and genetic differentiation 
in our target site suggest that relatively subtle alterations in habitat 
quality can strongly impact local population dynamics by decreas- 
ing individual numbers and disrupting gene flow at a very fine 
scale, independently of barriers that isolate habitats. As our study 
area was grazed by two different flocks (see Fig. 17 in Hochkirch et 
al. 2014), the observed demographic and genetic patterns may be 
the result of a difference in grazing pressure. 

This hypothesis was supported by the results of the analysis 
of vegetation productivity indices and NDVI in relationship with 
grazing treatment. Indeed, we found that three out of the four 
measured vegetation indices were significantly different between 
plots located inside and outside the fenced enclosure (Suppl. 
material 1: Fig. S6). The mean vegetation height and the cover- 
age in forbs were significantly higher within the fenced enclosure 
(Wilcoxon test p-value = 6.3e°’ and 0.0002, respectively) while 
the coverage in bare ground was significantly lower within the 


fenced enclosure (Wilcoxon test p-values = 0.0001). No difference 
was observed for the coverage in dry vegetation (Wilcoxon test 
p-values = 0.2302). The 2016-rescaled NDVI values were signifi- 
cantly correlated to the three productivity indices that responded 
to grazing. A positive correlation was found for the mean vegeta- 
tion height (Rho = 0.44; p-value = 6.4e°°) and the coverage in 
forbs (Rho = 0.33; p-value = 0.0008) while a negative correlation 
was found with the coverage in bare ground (Rho = -0.29; p-value 
= 0.004). Finally, the temporal analysis of the rescaled NDVI val- 
ues revealed much higher values for the vegetation plots located 
within than outside the enclosure in May 2016, after the fencing 
(Wilcoxon test p-value < 2.2e'°), while a slight opposite trend was 
found in May 2014, before the fencing (Wilcoxon test p-value = 
0.0563) (Suppl. material 1: Fig. S7). Within the area delineated by 
the fenced enclosure, rescaled NDVI values associated to vegeta- 
tion plots were also higher in 2016 than in 2014 (Wilcoxon test 
p-value < 2.2e'°). 

Altogether, this study suggests that P. rhodanica is sensitive to 
habitat quality and complements previous findings of a low disper- 
sal capability at the scale of the fragmented landscape. This may ex- 
plain why some subpopulations are no longer detected in the Crau 
Plain and imply that the few remaining ones may become extinct in 
the long-term as they are unlikely to be rescued through immigra- 
tion. This finding emphasizes the need for managing the P. rhodanica 
population at a local scale by considering the quality of the relict 
habitat patches, in addition to habitat fragmentation at a larger scale 
(i.e. sizes of and distances between Coussoul patches). Although 
this study did not identify clearly the processes driving this critically 
endangered species to extinction, the MAPI correlative approach 
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Fig. 6. Maps of A. density of grasshopper in number of individuals per hectare, B. rescale NDVI values and C. mean genetic differentia- 
tion between individuals resulting from the MAPI analysis. 
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Fig. 7. Spatial cross-correlograms between A. grasshopper density 
and NDVI, B. grasshopper density and the mean genetic differ- 
entiation between individuals (MAPI cell values) and, C. NDVI 
and the mean genetic differentiation between individuals. The x- 
intercept of the spline-correlogram is the estimate of the distance 
at which the correlation between variables is not different than 
expected by chance alone. Dotted lines represent the 95% confi- 
dence envelope based on 500 bootstrap resamples. 


helped us identify sheep grazing as a candidate landscape feature 
that may decrease grasshopper density and restrict gene flow within 
habitat patches. As our study was limited to a single sampling site, 
generalizing our results to the entire P. rhodanica population should 
be done with caution. Nonetheless, now that our indirect data-driv- 
en exploratory approach identified grazing pressure as a potential 
candidate driver of population decline, further work is needed in 
order to test for its population effects in a more direct way, draw firm 
conclusions and guide management actions. Above all, further fine 
monitoring of habitat quality (e.g. vegetation cover, structure and 
composition) in relation to direct measures of grazing pressure is 
critical. If the negative role of intense grazing is confirmed, imple- 
menting an adaptive management of pastoralism in the Crau Plain 
could help to sustain a higher number of reproductive grasshoppers 
and potential dispersers. 
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